A novel factor graph-based optimization technique for stereo correspondence estimation

Dense disparities among multiple views are essential for estimating the 3D architecture of a scene based on the geometrical relationship between the scene and the views or cameras. Scenes with larger extents of homogeneous textures, differing scene illumination among the multiple views and with occluding objects affect the accuracy of the estimated disparities. Markov random fields based methods for disparity estimation address these limitations using spatial dependencies among the observations and among the disparity estimates. These methods, however, are limited by spatially fixed and smaller neighborhood systems or cliques. Recent learning-based methods generate rich set of stereo features for generating cost volume and estimating disparity. In this work, we present a new factor graph-based probabilistic graphical model for disparity estimation that allows a larger and a spatially variable neighborhood structure determined based on the local scene characteristics. Our algorithm improves the accuracy of disparity estimates in stereo image pairs with varying texture and illumination characteristics by enforcing spatial dependencies among scene characteristics as well as among disparity estimates. We evaluated our method using the Middlebury benchmark stereo datasets and the Middlebury evaluation dataset version 3.0 and compared its performance with recent state-of-the-art disparity estimation algorithms. Our factor graph-based algorithm provided disparity estimates with higher accuracy when compared to the recent non-learning- and learning-based disparity estimation algorithms. The factor graph formulation can be used for obtaining maximum a posteriori estimates from models or optimization problems with complex dependency structure among hidden variables. The strategies of using a priori distributions with shorter support and spatial dependencies were useful for reducing the computational cost and improving message convergence in the model. The factor-graph algorithm is also useful for other dense estimation problems such as optical flow estimation.

(VSW) 17 , using adaptive support weight (ASW) 18 , and using a cross-based support window 19 . Using edge-preserving filters such as a bilateral filter (BF) 18,20 and guided image filter (GIF) 21 for cost aggregation provides a significant improvement on the final disparity estimates over the initial disparity estimates. An optimal disparity map or disparity estimates are obtained from the aggregated cost volume using optimization procedures such as the winner-take-all (WTA) 22 or using global optimization procedures such as dynamic programming(DP) 23 , graph cut (GC) 24 , intrinsic curves 25 , and belief propagation 26

methods.
Recently, learning-based methods have been successful in estimating stereo disparity with higher accuracy. Among the recent learning-based disparity estimation techniques, FC-DCNN is a light-weight, fully convolutional and densely connected neural network 27 based on a densely connected convolution network architecture 28 and a CNN based patch-matching technique for estimating stereo disparity 29 . FC-DCNN is a fully convolutional siamese neural network with individual network branches with shared weights for simultaneously processing left and right stereo images. With 5 convolutional layers and 64 filters in each layer, the network generated features were used to estimate similarity between a reference patch centered at each pixel in the left image and candidate patches from the right image in the stereo pair. For each pixel location, a cosine similarity measure was estimated using the corresponding feature vectors from the left and right channels of the network. The cost volume was filtered along all three dimensions and an initial disparity map was obtained using a winner-takes-all strategy. Accuracy of the initial disparity estimate at each pixel was assessed using a left-right consistency check. Inconsistent locations were identified and removed from the disparity map. Disparity values at these inconsistent locations were estimated based on classification of these locations as either part of homogeneous scene background or textured foreground objects, and using additional morphological processing of the disparity map. When compared to the MC-CNN-ACRT 29 , GC-Net 30 , and PSM-Net 31 , FC-DCNN is a lightweight network with 0.37 M trainable parameters with higher post-processing computational cost.
A successful strategy for improving accuracy of disparity estimates is to utilize spatial dependencies of the scene characteristics as well as of the disparity estimates. Among the probabilistic inference formulations for estimating dense geometric correspondences 32 , Markov Random Fields (MRF) based approaches have been successful in modeling spatial geometrical dependencies 33 . In brief, unknown true disparities among multiple views of a scene are modeled as random variables on the pixel lattice (random field) and dependence among random variables is assumed to follow Markov property. Thus, MRF represents a joint distribution of the random field using conditional distributions of each of the random variables. In learning-based MRF models, parameters of the MRF potential functions are either learned separately from training data 34,35 or along with the unknown states of the random variables 36 .
One of the limitations of the MRF models is that the neighborhood system used for enforcing spatial dependencies needs to be maximal. Further, the chosen dependency structure is uniformly enforced for all the random variables. Therefore, pairwise cliques or 2 × 2 cliques are more commonly used in MRF models. While learning methods are available for optimizing MRF parameters and neighborhood structure, they are generally limited to specific tasks.
Establishing dense correspondences among geometrical coordinates of multiple views is an ill-posed problem due to scene occlusion. Occluded regions are locations that are visible only in one of the images in a stereo image pair. Therefore, it is not possible to directly estimate disparities in these locations. Other primary sources of errors in disparity estimation are presence of larger scene extents with smoother texture characteristics, differences in scene illumination with respect to the multiple views/observers, and presence of sharp depth or disparity boundaries.
In this paper, we present a new factor graph-based probabilistic graphical model (FGS algorithm) that provides more accurate estimates of stereo disparity by using spatial dependency of scene characteristics (image intensity) and of disparity estimates. Figure 1 shows a schematic representation of computational steps in the FGS algorithm for generating a posterior disparity estimates for each stereo image pair. From a reference (left) image in each stereo pair, locations of various scene elements (foreground objects, and background elements) were grouped into image segments. A sparse but highly confident set of disparity estimates within each segmented region were used to build a priori distribution of disparities within each segmented region. In a factor graph model, all pixel locations within an image segment were associated with the respective a priori disparity distribution. For each pixel location, a posteriori disparity was estimated by enforcing spatial dependency among disparity estimates using a larger but spatially variable neighborhood system. A larger neighborhood system improves disparity estimates in regions with smoother texture and occlusion. A spatially variable neighborhood system provides higher disparity contrast along the depth boundaries. A final disparity map was obtained by post-processing the a posterior disparity map generated by the factor graph.
The main contributions of this study are as follows: • We present a new factor graph-based disparity estimation algorithm. To the best of our knowledge, this is the first application of factor graphs for estimating stereo and multi-view disparity. • The probabilistic graphical model allows use of larger and spatially variable neighborhood structure for enforcing belief exchange and dependency among highly related neighboring pixel locations. While larger neighborhood systems are beneficial for improving disparity estimates, the computational cost of belief exchange among neighbors in a graphical model is prohibitively high. Our FGS framework allows use of larger neighborhood systems with reduced computational steps by retaining disparity-dependency relationship only among highly related neighboring locations. • The FGS algorithm generates disparity maps with higher contrast along depth boundaries. This is possible because of our use of spatially-variable dependency structures while retaining dependency only with highly www.nature.com/scientificreports/ related neighboring locations (e.g. by discarding dependencies with pixels across an object or depth boundary). • We present strategies for reducing computational cost and for accelerating convergence of belief message exchanges in the factor graph namely 1. by reducing computational steps required for calculating belief measures in the factor graph by using binary potential functions; 2. by reducing number of variables as well as their limits of summation for marginalizing joint distributions during belief propagation in the factor graph; and 3. by generating a priori distributions with shorter support. • The factor graph structure can be used for determining a posteriori estimates in optimization problems with complex dependency structure among hidden variables.
We demonstrate the performance of the proposed probabilistic factor graph model by conducting extensive experiments using the Middlebury benchmark stereo datasets [37][38][39][40] and compare its performance with other state-of-the-art disparity estimation algorithms using Middlebury evaluation dataset version 3.0 5 . The remainder of the paper is structured as follows. In Section "Probabilistic factor graph model for disparity estimation", we present a detailed description of the new probabilistic factor graph-based stereo disparity estimation (FGS) algorithm. We present our experimental results in Section "Experimental results and discussion" and conclude this work in Section "Conclusions".

Probabilistic factor graph model for disparity estimation
For notational convenience, a linear index i ∈ {1, . . . , MN} was used to identity each of the pixels in images of size M × N pixels. Rectified stereo image pairs were corrected for any illumination differences using a homomorphic filter 41 . In brief, each image I(i) is modeled as an interaction between an illumination component l(i) and a reflectance component r(i) as I(i) = l(i) r(i) . Assuming that the illumination varies gradually over the imaging area, the illumination variation is subtracted from each image in the log domain using a high-pass filter. For each of the rectified and illumination corrected stereo image pairs, disparities between the left and right images were calculated by using the left image as reference.
Graph structure and message passing for approximate inference. A list of symbols used in the remainder of the paper is presented in Table 1. Figure   www.nature.com/scientificreports/ E = {1, . . . , MN} provide prior beliefs or evidences in assigning possible disparity labels to pixel locations. Dependency factor nodes S = {1, . . . , MN} are used to model spatial dependencies among the disparity labels assigned to the neighboring pixels. A random variable d i ∈ D is assigned to each variable node i ∈ V to represent the disparity label assigned to the ith pixel location. Each jth evidence factor node in E is connected one-to-one with the corresponding ith variable node in V to incorporate prior belief or evidence in determining the disparity labels d i . To represent the influence of each pixel location on its neighboring pixels, each variable node i ∈ V is connected to one or more dependency factor nodes {k ∈ S} using localized intensity characteristics in the reference image as described in Section "Determining variable nodes associated with each dependency factor node".
Let, ne(i) represent the set of neighboring nodes connected with any given node i; ne(i) \ j represent a set of all neighboring nodes of i excluding node j, where \ represents the set subtraction operation; D j = d i , ∀i ∈ V : i ∈ ne j ∈ F be a collection of random variables associated with any factor node j ∈ F ; A collection of random variables associated with any factor node j ∈ F An evidence potential function associated with factor node j ψ k , k ∈ S A dependency potential function associated with factor node k D Maximum number of disparity levels in a stereo pair Q Maximum number of variables nodes connected to a dependency factor node www.nature.com/scientificreports/ and let, D j \ i ⊆ D be a collection of random variables in D j except d i . An evidence potential function ψ j associated with each factor node j ∈ E is defined as a function of random variables D j of its neighboring nodes ψ j D j = ψ j d j . Similarly, the potential function ψ k (D k ) at the kth dependency factor node k ∈ S is defined as a function of the random variables D k associated with its neighboring nodes. Therefore, the factor graph represents the joint distribution of disparity labels assigned to each of the pixel locations as where, Z is the partitioning function. Probability of assigning various disparity labels to each pixel i can be obtained by marginalizing Eq. (1) with respect to D \ i.
This provides a sum-product formulation 42,43 for determining likely disparity labels for each of the pixel locations i based on a priori disparity information and spatial dependency characteristics of disparities in a stereo image pair.
For an approximate and efficient computation of marginal beliefs or probabilities in Eq. (2) using loopy belief propagation, local information available in each node is shared with neighboring nodes as variable-to-dependency factor messages µ i∈V→f ∈S and factor-to-variable messages µ f ∈F→i∈V until convergence 44,45 . Each outgoing message from a node is defined as a function of incoming messages at the given node as follows 44 .
Probabilistic model. For approximate inference on disparity label assignment, message structures for loopy belief propagation in Eqs. (3) and (4) and relevant potential functions were defined based on the posterior probability of disparity label assignment, where, f j D j is a function of the state of all variable nodes neighboring the factor node j. In our proposed method, the most relevant and compact spatial dependencies are defined independently at each pixel i using local image characteristics as described in Section "Determining variable nodes associated with each dependency factor node". Therefore, the joint likelihood term can be simplified as and the posterior probability updated as It can be observed that the posterior probability of d i without conditioning on a dependency factor node k ∈ ne(i) resembles the structure of the variable-to-dependency factor node message µ i→k in Eq. (3). This further suggests that individual likelihood terms p f j D j | d i form the dependency factor node-to-variable node message structure in Eq. (4).
Considering the individual likelihood terms in Eq. (6), www.nature.com/scientificreports/ and assuming that random variables D j associated with a dependency factor node j ∈ S are independent as in Moon & Gunther 46 , In the absence of any loops between any two variable nodes k and i (i.e. when there is at most one common dependency factor node s ∈ S between any two variable nodes k and i), the conditional probability p(d k | d i ) can be interpreted as the posterior probability of d k conditioned on all dependency factor functions associated with the kth variable node, i.e.
, the individual likelihood expression in Eq. (7) resembles the dependency factor node-to-variable node message structure µ j→i in Eq. (4).

Message passing implementation.
Based on the message-passing structures in our factor graph model, the variable-to-dependency factor node message µ i∈V→f ∈S in Eq. (3) was approximated as the posterior probability of d i conditioned on (satisfying) all of its neighboring factor nodes except the factor node f to which the message is sent as in Eq. (6). Similarly, the factor node-to-variable node messages µ f ∈F→i∈V in Eq. (4) was approximated as the likelihood of satisfying spatial dependency among the random variable states of all variables nodes associated with the factor node f except d i .
It can be observed that, for evidence factor nodes f ∈ E , the factor-to-variable messages in Eq. (4) simplifies as µ f ∈E→i∈V = ψ f (d i ) . This supplies fixed prior information about the state d i of the ith variable node by restricting messages from variable nodes to only dependency factor nodes as in Eq. (3). We estimated the a priori distribution p(d i ) from a disparity cost volume containing the cost of assigning all possible disparity labels at each pixel location. Details of building the disparity cost volume is presented in Section 2.4. Therefore, the approximate inference obtained using our factor graph model provides updated disparities D (an optimal surface within the cost volume) based on their posterior probabilities.
The potential function ψ f D f in Eq. (4) and its probabilistic representation p f j D j in Eq. (7) at the spatial dependency factor node j was assigned a value of 1.0 when the states D j \ i of the neighboring nodes are same as that of d i ; and was assigned a value of 0.0 when the states D j \ i are different from that of d i . This enforces spatial dependencies among the neighboring pixels in the final disparity map D . In addition, this further reduces the number of marginalization operations in Eqs. (4) and (7).
For each stereo image pairs, each evidence factor nodes e ∈ E were initialized with a priori probabilities p(d i ) of the variable node i ∈ ne(e) as the evidence factor node-to-variable node messages. The initial message from other variable and dependency factor nodes were set to be a uniform probability vector representing equally likely states. After initial message passing from the evidence factor nodes E , message exchange continues among all graph nodes V ∪ F until message convergence. We utilized an L 2 measure of for assessing message convergence at message passing iteration t + 1.

Disparity cost volume and a priori disparity distribution p(d i ).
For any given stereo image pair, let, C(i, d i ) be a cost volume representing the cost of assigning a disparity label d i to location i. Thus the final disparity map for a given stereo image pair will be an optimal surface within the cost volume C(i, d i ) . Figure 3 shows a schematic representation of algorithmic steps used for disparity cost volume calculation.
For computationally efficient disparity estimation, the reference image was segmented using an unsupervised texture segmentation method 47 . In brief, sharp image segment boundaries were derived using a Gabor filter bank and image boundaries were aggregated using k-means clustering to generate an image segmentation map. Within each of the segmented regions, highly confident disparity estimates at several candidate locations were obtained using an eigen-based feature matching method 48 . Using the zonal/regional disparity distributions, disparity cost at each of the location i was calculated using a normalized cross-correlation measure in the frequency domain. Within each segmented region, only disparity labels ranging from [µ d − σ d , µ d + σ d ] were considered based on the distribution of the disparity labels within the segmented region as shown in Fig. 3. This results in a sparse disparity cost volume C(i, d i ) and thus facilitates a faster inference due to reduced marginalization limits in Eq. (2). A detailed description of the initial cost volume computation as part of a hybrid of cross-correlation and scene segmentation (HCS) algorithm is available elsewhere 49 .
A priori probability of assigning a disparity label d i to the ith pixel location in the reference image was estimated from the cost volume C(i, d i ) as www.nature.com/scientificreports/ Determining variable nodes associated with each dependency factor node. We utilized edgepreserving filter kernels, namely the guided image filters (GIF) 50 and bilateral filters (BF) 51 to determine neighboring variable nodes with the highest influence on the true state of disparity at each of the ith pixel locations. Neighborhood dependency information from these kernels were used to define a non-symmetric, irregular, and higher-order neighborhood system based on the fact that objects at various depths in the scene may exhibit a disparity boundary along the object boundary. Let x, y represent the 2D coordinates of the ith pixel. Highly dependent neighbors of location (x, y) were selected using an αth percentile cut-off of the kernel coefficients. Coordinates of these highly dependent neighbors were used to identify the neighboring variable nodes of each dependency factor node. In brief, guided image filtering (GIF) is an edge-preserving smoothing algorithm. At each pixel location (x, y) in the reference image r l , a guided filter kernel W xy i, j with smoothness parameter ǫ was estimated as where ω xy is the window size used for estimating local illumination characteristics namely the mean illumination µ xy and standard deviation σ at location x, y .
The Bilateral filter (BF) 51 is an edge-preserving non-linear Gaussian filter with coefficients defined as a function of spatial and intensity similarities estimated respectively using a localized domain kernel and a range kernel. Bilateral filter kernel coefficients at location x, y is given as where the first exponential term represents a domain kernel as a function of pixel distance with respect to the center pixel x, y and the second term represents a range kernel as a function of regional image intensity with respect to that of the center pixel x, y . Parameters σ d and σ r controls the extent of influence neighboring pixels have on the domain and range filters respectively. Disparity estimates. Upon message passing convergence, an approximate estimate of the posterior probability of assigning a disparity label d i is available in each variable node as given in Eq. (5). A maximum a posterior disparity estimate was determined at each pixel location i based on the disparity label d i with the maximum posterior probability at the respective pixel i.
Post-processing the disparity maps. Occluded pixels were identified based on a lack of consistency or agreement between the disparity maps estimated using the left-right order vs right-left order of each stereo pair 52 . For each occluded pixel, the disparity estimate from its nearest non-occluded pixel within the same scanline (row) was assigned. Further, a weighted median filter was used to minimize spurious disparity assignments in the occluded region 53 . www.nature.com/scientificreports/

Experimental results and discussion
Algorithmic steps of the proposed factor graph method are presented in Algorithm 1. We evaluated the proposed method using stereo images from the Middlebury benchmark stereo datasets [37][38][39] , and 40 . We present a detailed evaluation of the proposed FGS algorithm followed by a comprehensive comparison of its performance with other state-of-the-art disparity estimation algorithms using Middlebury evaluation dataset version 3.0 5 .

FGS parameters and FGS implementation. The FGS algorithm was implemented in MATLAB 2018b
and its performance was evaluated using an Intel(R) Xeon(R) workstation with E3-1271 v3, 3.6 GHz processor. For illumination correction, high pass filters of size 21 × 21 pixels were obtained from a low-pass averaging filter. For unsupervised texture segmentation of the reference image, Gabor filters were designed 47 with filter orientations between 0 • − 150 • degrees in steps of 30 • , and wavelength starting from 2.83 up to the magnitude of hypotenuse of the input image. These filter parameters were previously identified to provide localized frequency features (roughly orthogonal subset) and spatial orientation features from input images 47 . K-means clustering algorithm was initialized with K = 15 with 5 replicates and ran for a maximum of 500 iterations. Based on our experiments, we observed that a higher number of clusters K are beneficial in reducing the computational complexity and for faster convergence of the FGS algorithm when there are several foreground objects at multiple depths in the scene. This is because of the fact that disparity values of locations within each distinct scene element are likely to be similar. Therefore, more accurate grouping of these scence elements will likely provide a priori disparity distributions with smaller standard deviations and shorter support leading to faster convergence of www.nature.com/scientificreports/ the FGS algorithm. While computational cost may be higher with non-optimal choice of clustering parameters the accuracy of the FGS disparity estimates is not affected. We evaluated various image features to identify robust features necessary for building a priori disparity distributions. Among the minimum eigenvalue algorithm, speeded up robust features algorithm (SURF), Harris corner and edge detector algorithm, and accelerated segment test algorithm (FAST), we observed that the eigen-based feature matching method 48 performed well with various types of benchmark images. For cost volume calculations, an optimal template size of 3 × 3 pixels were used 49 . For identifying variable nodes connected to each dependency factor node, bilateral filter kernel of size 7 × 7 pixels, domain kernel parameter of σ d = 3 , range kernel parameter of σ r = 0.1 and coefficient percentile cut-off of α = 97 were used. We observed that the smallest kernel sizes along with a higher percentile cut-off α identified fewer but highly dependent neighboring nodes. Therefore, the computational cost of message passing was significantly reduced with fewer but highly reliable neighboring variable nodes connected to each dependency factor node.
Computational complexity. In the FGS algorithm, message exchanges between the variable nodes V and the spatial dependency factor nodes S with |V| = |S| = MN are computationally demanding. Let, D be the cardinality of sample space of the random variable d i associated with the ith variable node i.e. maximum number of disparity levels in a stereo pair. Let, Q be the maximum number of variable nodes connected to a dependency factor node. Preparing a message for sending from a factor node f to a variable node i, µ f ∈F→i∈V in Eq. (4), requires marginalizing for the variable d i from the product form of joint distribution of disparity labels. The joint distribution is comprised of product of a potential function with (Q − 1) variable-to-factor messages. Factor node message related to d i is obtained by marginalizing a maximum of (Q − 1) neighboring variables D f \ d i . For each variable being marginalized, there are Q product operations i.e. ( Q − 1) + 1 . Because there are at most D disparity levels, marginalizing each variable requires D summation terms. Therefore, there are (D Q−1 Q) computations for each level of the disparity variable d i . With D levels of d i , there are D Q Q calculations for preparing a dependency factor-to-variable node message. In contrast, no calculations are required for preparing a message for sending from an evidence factor node to its corresponding variable node. Therefore, the FGS algorithm has a polynomial time complexity of O MN D Q Q per iteration including the computations required for preparing messages from each of the MN dependency factor nodes. It can be observed that far fewer number of calculations are required for preparing a message for sending from a variable node µ f ∈F→i∈V in Eq. (4) with a complexity of O(MN D (Q − 1)) and therefore does not affect the overall complexity of the FGS algorithm.
As described earlier, the following key features were employed to further reduce the computational requirements of the FGS algorithm.
1. (D Q Q) number of computations are eliminated when the state of the ith variable node d i is different from those of its neighboring nodes by using a binary potential function ψ f D f in the factor node-to-variable node messages µ f ∈F→i∈V in Eq. (4). Interaction among neighboring nodes can be further increased with appropriate changes to the potential function. 2. The number of disparity levels D in each of the (Q − 1) marginalization operations required for preparing factor node messages are significantly lowered by using a sparse cost volume for estimating a priori disparity in Eq. (9). 3. The number of factor node messages were reduced by connecting only highly correlated or dependent variable nodes (with a high bilateral filter coefficient cutoff) to each of the spatial dependency factor nodes.
Performance metrics. For performance evaluation, we used the common performance metrics available for assessing the accuracy of the estimated disparity maps namely, the disparity error maps, peak signal-tonoise ratio (PSNR), and average absolute error (Avg. err in pixels). Disparity error maps were computed as location-wise difference between the estimated disparity d (x, y) and its ground-truth d(x, y) as d (x, y) − d(x, y) . PSNR provides a measure of similarity between an estimated disparity map d (x, y) of size M × N pixels and the ground-truth disparity map d(x, y) as follows.
A thresholded average disparity error metric with a disparity threshold of T pixels was defined as Average disparity errors were assessed at two disparity threshold levels of T = 2 pixels (Bad2.0) and T = 0 pixels (Avg. err).
In the majority of the stereo pairs tested, the FGS algorithm converged between 25 and 30 iterations based on the convergence measure in Eq. (8).
Filter selection for identifying FGS variable node neighbors. The edge-preserving filters (Sect. 'Determining variable nodes associated with each dependency factor node') with the highest accuracy and computational speed were selected for identifying neighboring variable nodes for each of the dependency (12) PSNR = 10 log 10 www.nature.com/scientificreports/ factor nodes in the FGS algorithm. Figure 4c and e show the disparity maps for the "Teddy" stereo pair (ground truth disparity shown in Fig. 4a) estimated using guided image filters and bilateral filters respectively based on an initial disparity shown in Fig. 4b. Disparity error maps for the guided image filter and bilateral filter are shown in Fig. 4d and f respectively. Quantitative performance measures of the edge-preserving filters are presented in Table 2. The accuracy of the FGS algorithm with guided filters (based on the Avg. err, PSNR, and Bad2.0 metrics) were slightly better than the bilateral filters. Because the guided filters resulted in a larger number of neighboring variable nodes, the run time of the FGS algorithm, however, was higher with the guided filters when compared to the FGS algorithm with bilateral filters. Therefore, for further evaluation of the FGS algorithm, bilateral filtering was chosen as the optimal choice for identifying neighboring variable nodes in the FGS factor graph. Assessment results without post-processing. For the selected stereo pairs of scenes with differing textures and scene illumination, Fig. 5 shows disparity maps estimated by the HCS algorithm without cost aggregation and post-processing, disparity maps estimated by the FGS algorithm without any post-processing, and errors in the disparity maps estimated by the FGS algorithm. A summary of the assessment metrics without post-processing the disparity maps is presented in Table 3.
Without post-processing, the FGS algorithm provided significantly lower average error (Avg. err) and percentage of bad pixels (Bad2.0) when compared to the initial disparity estimates. Higher PSNR values indicate a higher degree of similarity between the FGS estimated disparity maps and the ground truth disparity maps. Visual  www.nature.com/scientificreports/ inspection of the FGS disparity maps in Fig. 5 revealed that the FGS algorithm is able to estimate more accurate disparities in occluded regions and generate disparity maps with well-defined boundaries. These observations validate our choice of using non-symmetric / irregular, variable and higher-order neighborhood dependency relationship among neighboring pixel locations in the FGS algorithm.
Assessment results with post-processing. For selected stereo pairs with varying texture and illumination characteristics, Fig. 6 shows the disparity maps generated after post-processing the maximum a posteriori dispar-  www.nature.com/scientificreports/ ity estimates (shown in Fig. 5) and their corresponding disparity error maps. Table 4 presents a summary of assessment metrics after post-processing the disparity estimates. After post-processing, inconsistencies in stereo disparity estimates were corrected and disparity maps with sharper depth boundaries were generated. The average error (Avg. err) and percentage of bad pixels (Bad2.0) further decreased and the PSNR value increased after post-processing indicating higher similarity with the ground truth disparity estimates. Though the FGS disparity estimates without post-processing demonstrated higher accuracy in occluded and homogeneous regions, postprocessing using a weighted median filter further reduced inconsistent disparity estimates and improved the overall quality of the disparity maps.

Performance of the FGS algorithm vs state-of-the-art algorithms.
To the best of our knowledge, the proposed FGS algorithm is the first disparity estimation technique based on factor graphs. The performance of the non-learning-based FGS algorithm was compared with recent state-of-the-art learning-based and nonlearning-based disparity estimation algorithms. All disparity estimation algorithms were evaluated using stereo pairs from the current Middlebury evaluation dataset version 3.0. In addition to the performance metrics of Avg. err, PSNR and Bad 2.0 %, a weighted average measure was calculated for each performance metric based on the level of difficulty of estimating the disparity map for each stereo pair in the Middlebury evaluation dataset ver-  www.nature.com/scientificreports/ sion 3.0. A summary of all the assessment metrics and a weighted performance measure for the FGS algorithm is presented in Table 5.

FGS vs non-learning based disparity estimation methods.
Performance of the FGS algorithm was compared with 12 recently developed non-learning-based disparity estimation procedures including 7 local methods, 3 global methods, and one fusion method using both local and global approaches for disparity estimation. The local methods were: weighted adaptive cross-region-based guided image filtering method (ACR-GIF-OW) 54 , real-time stereo matching algorithm with FPGA architecture (MANE) 55 , adaptive support-weight approach in pyramid structure (DAWA-F) 56 , encoding-based approaches PPEP-GF 57 , absolute difference (AD) and census transform-based stereo matching with guided image filtering (ADSR-GIF) 58 , the sum of absolute difference (SAD) based stereo matching aggregated with adaptive weighted bilateral filter (SM-AWP) 59 , statistical maximum a posteriori estimation of MRF disparity labels (SRM) 60 . The global disparity estimation procedures chosen for comparison were: binocular narrow-baseline stereo matching procedure using a max-tree data structure (MTS) 61 and its improvement (MTS-2) 62 , and an accelerated multi-block matching (MBM) algorithm on GPU 22 . FASW approach 63 uses both local and global strategies and is based on a census transform with adaptive support weight. Avg. err metric and a weighted average of Avg. err for the FGS algorithm and the non-learning-based disparity estimation procedures are presented in Table 6. Image weight given in Table 6 represents the level of difficulty in estimating the disparity from a given stereo pair. In general, the proposed FGS algorithm provided higher accuracy for all the stereo pairs comparable to that of other non-learning-based methods. The FGS method provided the lowest weighted average of Avg. err of 6.45 pixels among all the non-learning-based methods. Further, the FGS method provided the lowest estimation error (Avg. err) for 3 out of the 10 difficult stereo pairs (image weight = 8) and for 4 out of the 5 moderately difficult stereo pairs (image weight = 4). Among the remaining 7 difficult stereo pairs, HCS method provided lower error for 1 pair, ACR-GIF-OW method for 1 pair, SRM method for 2 pairs, DAWA-F method for 1 pair, and FASW for 2 pairs. The DAWA-F method provided the lowest estimation error for the remaining 1 moderately difficult stereo pairs. FGS vs learning-based disparity estimation methods. Performance of the FGS algorithm was also compared with the following recently developed learning-based disparity procedures: a method based on a fusion of convolutional neural networks (CNN) and conditional random fields (LBPS) 64 , a fully convolutionaldensely connected neural network (FC-DCNN) 27 , a deep-learning assisted method to produce initial estimate with Semi-Global Block Matching method (SGBMP) 65 , a deep learning-based self-guided cost aggregation method (DSGCA) 66 , a stereo matching algorithm with a pretrained network and global energy minimization SIGMRF 67 , a multi-dimensional convolutional neural network (MSMD-ROB) 68 and a CNN-based network using ResNet (CBMBNet) 69 .
Avg. err metric for the FGS algorithm and the learning-based disparity estimation procedures are presented in Table 7. Among the recent learning-based disparity estimation procedures, the LBPS method provided the least weighted Avg. err of 5.51 for all stereo pairs, provided the least Avg. error for 4 out of 10 difficult stereo pairs www.nature.com/scientificreports/ and the least Avg. err for 4 out of 5 moderately difficult stereo pairs. The FGS method provided the second-lowest weighted average of Avg. err of 6.45 for all stereo pairs, provided the least Avg. err for 4 out of 10 difficult stereo pairs and the least Avg. error for 1 out of 5 moderately difficult stereo pairs. The CBMBNet provided the third lowest weighted Avg. err of 6.65 and the least Avg. err for 2 out of 10 difficult stereo pairs.

Conclusions
We have presented a new probabilistic factor-graph-based disparity estimation algorithm that improves the accuracy of disparity estimates in stereo image pairs with varying texture and illumination characteristics by enforcing spatial dependencies among scene characteristics as well as among disparity estimates. In contrast to  Table 7. Performance (Avg. err) of the FGS algorithm vs state-of-the-art learning-based disparity estimation algorithms for stereo pairs in the Middlebury evaluation dataset version 3.0. The weight of a stereo pair represents the level of difficulty in estimating disparity map from the respective stereo pair. The algorithm with the highest performance is highlighted in bold. www.nature.com/scientificreports/ MRF models, our factor graph formulation allows a larger as well as a spatially variable neighborhood system dependent only on the local scene characteristics. Our factor graph formulation can be used for obtaining maximum a posteriori estimates from models or optimization problems with complex dependency structure among hidden variables. The strategies of using a priori distributions with shorter support and spatial dependencies are useful for improving the computational speed of message convergence in factor graph-based inference problems. We rigorously evaluated the performance of the new factor-graph-based disparity estimation algorithm using Middlebury benchmark stereo datasets [37][38][39] , and 40 . Our experimental results indicate that the factor-graph algorithm provides disparity estimates with higher accuracy when compared to recent non-learning as well as learning-based disparity estimation algorithms using Middlebury evaluation dataset version 3.0 5 . The factorgraph algorithm may also be useful for other dense estimation problems such as optical flow estimation.